ProjektOpVar
Mateusz Strojek
2024-04-21
Modelowanie ryzyka kredytowego i operacyjnego - projekt
Zadanie 1
Celem tego zadania jest wyznaczenie operacyjnego VaR i ES. Aby to zrobić, to wyznaczone zostanę częstości strat dla każdego roku dla każdej kategorii. Zdarzenia oznaczają kolejno:
- Zakłócenia w działalności i awarie systemów (Buss_distr)
- Klienci, produkty i praktyki biznesowe (Com_Ban) - tutaj czyste domysły z #### prezentacji
- Uszkodzenia mienia materialnego (Damage)
- Praktyki zatrudnienia i bezpieczeństwo w miejscu pracy (Empl_Pract)
- Wykonanie, dostarczenie i zarządzanie procesem (Execut_Del)
- Oszustwa zewnętrzne (External Fr)
- Oszustwa wewnętrzne (Internal Fr)
Następnie, za pomocą testów zgodności, zostanie dopasowany najlepszy rozkład do częstości strat oraz do ich dotkliwości. Zostanie wygenerowanych 20000 scenariuszy i zostanie wyznaczony symulacyjny rozkład strat. Procedura wyznaxczania OpVaR i OpES zostanie powtórzona 200 razy. Został mi przydzielony plik22.csv.
Wyznaczenie częstości strat
W tej sekcji zostały wyznaczone częstości strat, pokazano jak kształtują się one w tabelce w określonym roku, przedstawiono na wykresach oraz obliczono średnią, wariancję oraz średnią dzieloną przez wariancję w celu późniejszego wyboru, jakimi rozkładami zostaną dopasowane poszczególne częstości.
library(dplyr)
library("tidyverse")
library(fitdistrplus)
library(actuar)
library(DT)
library(gridExtra)
library(BatchGetSymbols)
options(scipen = 0, digits = 4)
dane <- read.csv("plik22.csv")
czestosc <- dane %>%
group_by(linia) %>%
count(rok)
wartosci <- unique(dane$linia)
ramka_danych <- data.frame("Year"=unique(dane$rok))
for(i in 1:length(wartosci)){
x <- czestosc[czestosc$linia==wartosci[i],3]
ramka_danych <- cbind(ramka_danych,x)
}
colnames(ramka_danych) <- c("Year",wartosci)
datatable(ramka_danych)Histogramy częstości dla poszczególnej linii
Buss_Distr
Com_Ban
Damage
Empl_Pract
Execut_Del
External_Fr
Internal_Fr
Tabelka ze średnimi, wariancjami i ich ilorazem:
Można zauważyć, że zdarzeniem o najmniejszej częstotliwości jest zdarzenie “Buss_distr”. Zdarza się ono ok. 5 razy w roku. Zdarzenie “Com_Ban” występuje ono ok. 18 razy w roku, przy czym warto podkreślić, że wartości te wachają dość mocno, gdyż współczynnik zmienności wynosi ok. 33%. Zdarzenie “Damage” występuje średnio 8 razy rocznie. Drugim najczęstszym zdarzeniem jest “Empl_Pract”, ponieważ występuje ono średnio 131 razy rocznie. Pomimo tego, że wariancja wynosi aż 308.6, to wsp. zmienności wynosi ok. 13%, zatem można stwierdzić że częstości tego zdarzenia nie wahają się zbyt mocno. Zdarzenie “Execut_Del” występuje średnio 20 razy rocznie, a “External Fr.” odpowiadające oszustwom zewnętrznym 10 razy. Najczęstszym zdarzeniem okazały się oszustwa wewnętrzne. Rocznie można było zaobserwować je ok. 161 razy. Współczynnik zmienności wynosi ok. 13%, więc ponownie jest to zdarzenie, którego częstość nie waha się mocno.
Wnioskując po wynikach ilorazu średniej przez wariancję, postanowiono, że:
- “Buss_distr” - rozkład dwumianowy (wynik mniejszy od 1)
- “Com_Ban” - rozkład ujemny dwumianowy (wynik większy od 1)
- “Damage” - rozkład dwumianowy (wynik większy od 1)
- “Empl_pract” - rozkład ujemny dwumianowy (wynik mniejszy od 1)
- “Execut_del” - rozkład ujemny dwumianowy (wynik mniejszy od 1)
- “External_Fr” - rozkład dwumianowy (wynik większy od 1)
- “Internal_Fr” - rozkład ujemny dwumianowy (wynik mniejszy od 1)
W żadnym przypadku nie wybrano rozkładu Poissona, gdyż wyniki wariancji i wyniki średniej nie były podobne same w żadnym przypadku.
Dopasowanie częstości zdarzeń do poszczególnych rozkładów:
Parametr size przy dopasowaniu rozkładu dwumianowego został wyznaczony z układu równań:
Średnia = np
Wariancja = np(1-p)
dwu1 <- fitdist(ramka_danych[,2], distr="binom",
start = list(prob=0.5),
fix.arg = list(size=13))
dwu2 <- fitdist(ramka_danych[,4], distr="binom",
start = list(prob=0.5),
fix.arg = list(size=15))
dwu3 <- fitdist(ramka_danych[,7], distr="binom",
start = list(prob=0.5),
fix.arg = list(size=27))
uj1 <- fitdist(ramka_danych[,3], method="mme",distr="nbinom")
uj2 <- fitdist(ramka_danych[,5], method="mme",distr="nbinom")
uj3 <- fitdist(ramka_danych[,6], method="mme",distr="nbinom")
uj4 <- fitdist(ramka_danych[,8], method="mme",distr="nbinom")
ramkapvalue <-
data.frame(c(gofstat(dwu1)$chisqpvalue,
gofstat(dwu2)$chisqpvalue,
gofstat(dwu3)$chisqpvalue,
gofstat(uj1)$chisqpvalue,
gofstat(uj2)$chisqpvalue,
gofstat(uj3)$chisqpvalue,
gofstat(uj4)$chisqpvalue))
colnames(ramkapvalue) <- c("P.value")
rownames(ramkapvalue) <- c("Buss_Distr", "Damage", "External_Fr",
"Com_ban", "Empl_pract", "Execut_Del",
"Internal_Fr")
datatable(ramkapvalue)Można zauważyć, że p-value dla testu chi-kwadrat wynosi wszędzie powyżej 0.099. Jednak warto zwrócić uwagę, że wyniki te mogą być mylące. Nie został spełniony warunek oczekiwanej liczebności poszczególnych komórek tablicy korelacyjnej (powyżej 5).
Wybór rozkładów dla dotkliwości strat i jego dopasowanie
Sprawdzono, czy dane pasują do rozkładu Weibulla, logarytmiczno-normalnego, wykładniczego, wykładniczego, Pareto oraz rozkładu gamma. Następnie, przeprowadzono testy Kołmogorowa-Smirnowa, Andersona Darlinga oraz chi-kwadrat. Test Kołmogorowa najlepiej wychwytuje różnice w okolicach średniej/mediany, ale nie jest dokładny w ogonach rozkładu. Z drugiej strony, test Andersona Darlinga jest bardziej czuły na różnice w ogonach.
W tabelce znajdują się wyniki testów dla danego rozkładu. Każda strona odpowiada jednej linii.
############################################################
Buss_Distr <- dane[dane$linia==wartosci[1],3]
Com_Ban <- dane[dane$linia==wartosci[2],3]
Damage <- dane[dane$linia==wartosci[3],3]
Empl_Pract <- dane[dane$linia==wartosci[4],3]
Execut_Del <- dane[dane$linia==wartosci[5],3]
External_Fr <- dane[dane$linia==wartosci[6],3]
Internal_Fr <- dane[dane$linia==wartosci[7],3]
lista <- list(Buss_Distr,
Com_Ban,
Damage,
Empl_Pract,
Execut_Del,
External_Fr,
Internal_Fr)
pvalue <- list()
for(i in 1:length(lista)){
r1 <- fitdist(lista[[i]], "weibull")
r2 <- fitdist(lista[[i]], method="mme", "lnorm")
r3 <- fitdist(lista[[i]], method="mme", "exp")
r4 <- fitdist(lista[[i]], method="mme", "gamma")
r5 <- fitdist(lista[[i]], distr ="pareto", start = list(shape = 1, scale = 500))
test <- gofstat(list(r1,r2,r3,r4,r5),
fitnames = c("weibull", "lnorm", "exp", "gamma","pareto"))
chisqpvalue <- ifelse(test$chisqpvalue < 0.05, "rejected", "not rejected")
pvalue[[wartosci[i]]] <- data.frame("KSTest"=test$kstest,
"ADTest"=test$adtest,
"ChisqTest"=chisqpvalue,
"ChisqPvalue"=test$chisqpvalue)
}## Warning in cov2cor(varcovar): 'diag(.)' miał 0 lub wpisy NA; nieskończony
## rezultat jest wątpliwy## Warning in sqrt(diag(varcovar)): wyprodukowano wartości NaN## Warning in sqrt(1/diag(V)): wyprodukowano wartości NaN## Warning in cov2cor(varcovar): 'diag(.)' miał 0 lub wpisy NA; nieskończony
## rezultat jest wątpliwydatatable(t(as.data.frame(pvalue)), options = list(pageLength = 4))Wybrałem odpowiednie rozkłady:
“Buss_Distr” - weibull, gdyż wszystkie testy nie odrzuciły H0, że dane należą do rozkładu Weibulla i p-value z testu chi-kwadrat było największe,
“Com_Ban” - gamma, ponieważ test kołmogorowa nie odrzucił H0 oraz wartość p-value z testu chi-kwadrat było największe spośród rozkładów, które nie odrzuciły H0 w teście Kołmogorowa,
“Damage” - weibull, ponieważ wszystkie testy nie odrzuciły H0 oraz p-value testu chi-kwadrat było większe od tego dotyczącego rozkładu gamma,
“Empl_Pract” - weibull, ponieważ wartość p-value testu chi-kwadrat było największe. Warto zwrócić uwagę że wszystkie testy dla każdego rozkładu odrzuciły hipotezę główną, że dane należą do danego rozkładu,
“Execut_Del” - weibull, sposób wyboru został podjęty w sposób analogiczny do zdarzenia “Empl_Pract”,
“External_Fr” - weibull, ponieważ wartość p-value testu chi-kwadrat było większe niż dla rozkładu gamma. Porównano wartości z tym rozkładem, ponieważ w przypadku testu KS oraz chi-kwadrat, nie odrzucono hipotezy głównej.
“Internal_Fr” - lnorm, ponieważ test KS nie odrzucił hipotezy H0. We wszystkich innych testach zarówno dot. rozkładu lnorm, jak i innych, odrzucano H0.
param_buss <- fitdist(lista[[1]], "weibull")
param_com <- fitdist(lista[[2]], method="mme", "gamma")
param_dmg <- fitdist(lista[[3]], "weibull")
param_empl <- fitdist(lista[[4]], "weibull")
param_exe <- fitdist(lista[[5]], "weibull")
param_ext <- fitdist(lista[[6]], "weibull")
param_int <- fitdist(lista[[7]], method="mme", "lnorm")
wyk1 <- ggplot(data = data.frame(lista[[1]]), aes(x = lista[[1]])) +
geom_histogram(aes(y = ..density..), color = "black", fill = "white", bins = 14) +
geom_density(aes(color = "Empiryczny")) +
geom_function(aes(color = "Weibull"),
fun = dweibull,
args = list(param_buss$estimate[1],
param_buss$estimate[2])) +
labs(title = "Buss_distr", x = "Buss_distr") +
theme(legend.position = "top") +
scale_color_manual(values = c("blue", "red"), labels = c("Empiryczny","Weibull")) +
theme_minimal()
wyk2 <- ggplot(data = data.frame(lista[[2]]), aes(x = lista[[2]])) +
geom_histogram(aes(y = ..density..), color = "black", fill = "white", bins = 14) +
geom_density(aes(color = "Empiryczny")) +
geom_function(aes(color = "Gamma"),
fun = dgamma,
args = list(param_com$estimate[1],
param_com$estimate[2])) +
labs(title = "Com_ban", x = "Com_ban") +
theme(legend.position = "top") +
scale_color_manual(values = c("blue", "red"), labels = c("Empiryczny","Gamma")) +
theme_minimal()
wyk3 <- ggplot(data = data.frame(lista[[3]]), aes(x = lista[[3]])) +
geom_histogram(aes(y = ..density..), color = "black", fill = "white", bins = 14) +
geom_density(aes(color = "Empiryczny")) +
geom_function(aes(color = "Weibull"),
fun = dweibull,
args = list(param_dmg$estimate[1],
param_dmg$estimate[2])) +
labs(title = "Damage", x = "Damage") +
theme(legend.position = "top") +
scale_color_manual(values = c("blue", "red"), labels = c("Empiryczny","Weibull")) +
theme_minimal()
wyk4 <- ggplot(data = data.frame(lista[[4]]), aes(x = lista[[4]])) +
geom_histogram(aes(y = ..density..), color = "black", fill = "white", bins = 14) +
geom_density(aes(color = "Empiryczny")) +
geom_function(aes(color = "Weibull"),
fun = dweibull,
args = list(param_empl$estimate[1],
param_empl$estimate[2])) +
labs(title = "Empl_Pract", x= "Empl_Pract") +
theme(legend.position = "top") +
scale_color_manual(values = c("blue", "red"), labels = c("Empiryczny","Weibull")) +
theme_minimal()
wyk5 <- ggplot(data = data.frame(lista[[5]]), aes(x = lista[[5]])) +
geom_histogram(aes(y = ..density..), color = "black", fill = "white", bins = 14) +
geom_density(aes(color = "Empiryczny")) +
geom_function(aes(color = "Weibull"),
fun = dweibull,
args = list(param_exe$estimate[1],
param_exe$estimate[2])) +
labs(title = "Execut_Del", x="Execut_Del") +
theme(legend.position = "top") +
scale_color_manual(values = c("blue", "red"), labels = c("Empiryczny","Weibull")) +
theme_minimal()+
xlim(c(0, 200))
wyk6 <- ggplot(data = data.frame(lista[[6]]), aes(x = lista[[6]])) +
geom_histogram(aes(y = ..density..), color = "black", fill = "white", bins = 14) +
geom_density(aes(color = "Empiryczny")) +
geom_function(aes(color = "Weibull"),
fun = dweibull,
args = list(param_ext$estimate[1],
param_ext$estimate[2])) +
labs(title = "Extern_fr", x="Extern_fr") +
theme(legend.position = "top") +
scale_color_manual(values = c("blue", "red"), labels = c("Empiryczny","Weibull")) +
theme_minimal()
wyk7 <- ggplot(data = data.frame(lista[[7]]), aes(x = lista[[7]])) +
geom_histogram(aes(y = ..density..), color = "black", fill = "white", bins = 14) +
geom_density(aes(color = "Empiryczny")) +
geom_function(aes(color = "Log-norm"),
fun = dlnorm,
args = list(param_int$estimate[1],
param_int$estimate[2])) +
labs(title = "Intern_fr", x="Intern_fr") +
theme(legend.position = "top") +
scale_color_manual(values = c("blue", "red"), labels = c("Empiryczny","Log-norm")) +
theme_minimal()Wykresy rozkładów
Buss_Distr
Com_Ban
Damage
Empl_Pract
Execut_Del
External_Fr
Internal_Fr
Widać widocznie, że dla zdarzenia “Execut_Del”, rozkład Weibulla nie do końca oddaje rzeczywisty rozkład. Z drugiej strony, rozkład dopasowany do dotkliwości strat zdarzenia “Empl_Pract” nie różni się tak bardzo od faktycznego rozkładu. Warto dodać, że wszystkie testy odrzuciły H0 dla tego zdarzenia.
Przeprowadzenie symulacji
Przeprowadzono symulację 200 razy. Za każdym razem wygenerowano 20000 scenariuszy. Aby to zrobić, wygenerowano dla każdej linii 20 tysięcy wartości z określonych rozkładów i połączono w jedną ramkę danych. Następnie wygenerowano x różnych wartości z określonego rozkładu dla poszczególnej linii. Za x odpowiada odpowiednia wartość w ramce danych określającej częstości strat. Następnie, po zsumowaniu dotkliwości strat ze wszystkich linii dla każdego scenariusza, obliczono VaR i ES.
liczba_zdarzen <- 20000
liczba_powtorzen <- 200
VaR <- c()
ES <- c()
for(i in 1:liczba_powtorzen){
Buss_czest <- rbinom(liczba_zdarzen, size = 13, prob = dwu1$estimate)
Dmg_czest <- rbinom(liczba_zdarzen, size = 15, prob = dwu2$estimate)
Ext_czest <- rbinom(liczba_zdarzen, size = 27, prob = dwu3$estimate)
Com_czest <- rnbinom(liczba_zdarzen,
size=round(uj1$estimate[1],0),
mu=uj1$estimate[2])
Empl_czest <- rnbinom(liczba_zdarzen,
size=round(uj2$estimate[1],0),
mu=uj2$estimate[2])
Exec_czest <- rnbinom(liczba_zdarzen,
size=round(uj3$estimate[1],0),
mu=uj3$estimate[2])
Int_czest <- rnbinom(liczba_zdarzen,
size=round(uj4$estimate[1],0),
mu=uj4$estimate[2])
ramka_czestosc <- data.frame("Buss"=Buss_czest,
"Com"=Com_czest,
"Damage"=Dmg_czest,
"Empl."=Empl_czest,
"Exec."=Exec_czest,
"Ext."=Ext_czest,
"Int."=Int_czest)
straty_ramka <- data.frame()
for(i in 1:nrow(ramka_czestosc)){
strata_buss <- sum(rweibull(ramka_czestosc[i,1],
shape=param_buss$estimate[1],
scale=param_buss$estimate[2]))
strata_com <- sum(rgamma(ramka_czestosc[i,2],
shape=param_buss$estimate[1],
scale=param_buss$estimate[2]))
strata_dmg <- sum(rweibull(ramka_czestosc[i,3],
shape=param_dmg$estimate[1],
scale=param_dmg$estimate[2]))
strata_empl <- sum(rweibull(ramka_czestosc[i,4],
shape=param_empl$estimate[1],
scale=param_empl$estimate[2]))
strata_exe <- sum(rweibull(ramka_czestosc[i,5],
shape=param_exe$estimate[1],
scale=param_exe$estimate[2]))
strata_ext <- sum(rweibull(ramka_czestosc[i,6],
shape=param_ext$estimate[1],
scale=param_ext$estimate[2]))
strata_int <- sum(rlnorm(ramka_czestosc[i,7],
meanlog = param_int$estimate[1],
sdlog = param_int$estimate[2]))
straty_ramka <- rbind(straty_ramka,
c(strata_buss,
strata_com,
strata_dmg,
strata_empl,
strata_exe,
strata_ext,
strata_int)
)
}
colnames(straty_ramka) <- colnames(ramka_czestosc)
straty_ramka$Suma <- rowSums(straty_ramka)
wartosci <- sort(straty_ramka$Suma, decreasing = T)
kwantyl <- 0.999
VaR <- c(VaR, quantile(wartosci, kwantyl))
indeks_do_wycięcia <- length(wartosci) * (1-kwantyl)
wyc_wartosci <- wartosci[1:indeks_do_wycięcia]
ES <- c(ES, mean(wyc_wartosci))
}Wygląd ramki odpowiadającej za częstości strat:
datatable(ramka_czestosc)## Warning in instance$preRenderHook(instance): It seems your data is too
## big for client-side DataTables. You may consider server-side processing:
## https://rstudio.github.io/DT/server.htmlWygląd ramki danych odpowiadającej za dotkliwość strat. Każdy rząd odpowiada za jeden scenariusz.
datatable(straty_ramka,
options = list(
columnDefs = list(list(
targets = "_all",
render = JS("function(data, type, row, meta) {",
"if(type === 'display' && typeof data === 'number') {",
"return parseFloat(data).toFixed(2);",
"} else {",
"return data;",
"}",
"}")
))
))## Warning in instance$preRenderHook(instance): It seems your data is too
## big for client-side DataTables. You may consider server-side processing:
## https://rstudio.github.io/DT/server.htmlObliczony VaR i ES:
datatable(
data.frame("VaR"= c(mean(VaR), sd(VaR), sd(VaR)/mean(VaR)), "ES" = c(mean(ES),sd(ES), 100*sd(ES)/mean(ES)),
row.names = c("Średnia","Odchylenie st.", "Współczynnik zmienności (%)"))
)Współczynnik zmienności obliczonych wartości wynosi max 1%, zatem można uznać wyniki za stabilne.
Zadanie 2
Czym jest model Mertona?
Model Mertona to tzw. model strukturalny, który wykorzystuje informacje z bilansów firm do modelowania ich aktywów jako instrumentów finansowych. W modelu tym zakłada się, że wartość aktywów firmy zmienia się według geometrycznego ruchu Browna – modelu Blacka-Scholesa, gdzie parametry są stałe. Przyjmuje się również, że ryzyko niewypłacalności firmy materializuje się tylko w momencie, gdy dojdzie do zapadalności wyemitowanych przez nią obligacji zerokuponowych. Często nie jest możliwe bezpośrednie obserwowanie bieżącej wartości rynkowej aktywów firmy oraz ich zmienności. W takich przypadkach, teoria opcji dostarcza matematycznych narzędzi do ich estymacji. Używając modelu Blacka-Scholesa, można obliczyć wartość kapitału własnego firmy w danym momencie, bazując na przewidywanej wartości i zmienności aktywów. Model ten pozwala na analizę wpływu zmienności rynkowej, przyszłej wartości aktywów, zobowiązań oraz stopy wolnej od ryzyka na wartość kapitału własnego.
Dla opcji kupna można wycenić opcję w chwili t=0 robiąc to za pomocą wzoru Blacka-Scholesa:
\[ E_t = A_tN(d_1) - Le^{-r(T-t)}N(d_2), \]
gdzie
\[ d_1 = \frac{\ln(A_t/L) + (r + (\sigma^2)/2)(T-t)}{\sigma\sqrt{T}} \]
oraz
\[ d_2 = d_1 - \sigma\sqrt{T} \]
Teoria opcji pozwala wyznaczyć prawdopodobieństwo, że firma nie spłaci tego zadłużenia. Liczy się ono następująco:
\(N(-d_2) = 1 - N(d_2)\).
Niestety, nie ma wystarczająco dużo informacji, które pozwalają wyznaczyć niektóre zmienne we wzorze. Z pomocą przychodzi lemat ito, czyli:
\[ \left\{ \begin{align*} E_0 &= V_0N(d_1) - De^{-rT}N(d_2), \\ \sigma_{E}E_0 &= N(d_1)\sigma_{V}V_0 \end{align*} \right. \]
Parametry można łatwo oszacować korzystając z minimalizacji błędu.
Wyznaczenie PD
W zadaniu użyto danych firmy Asseco Poland S.A. Jest to firma notowana na Giełdzie Papierów Wartościowych, pochodząca z branży IT. Wartość zobowiązań firmy wynosi 976500000 zł, co można ujrzeć w bilansie. Informacje o liczbie wyemitowanych akcji można zaczerpnąć ze strony https://notowania.pb.pl/instrument/PLSOFTB00016/assecopol/informacje-spolka. Jest ich razem 83000303. Stopę wolną od ryzyka ustanowiono na poziomie 6,25% (oprocentowanie rocznych obligacji skarbowych z dnia 1 grudnia 2023r).
Pobrano informacje odnośnie cen akcji na przestrzeni 2023 roku ze strony yahoo, co przedstawia kod poniżej:
first_date <- as.Date("2023-01-01")
last_date <- as.Date("2023-12-31")
l.out <- BatchGetSymbols(tickers = c("ASOZY"),
first.date = first_date,
last.date = last_date,
do.cache = FALSE)Obliczenie rocznej zmienności zwrotów logarytmnicznych a następnie ustalenie wartości kapitału. Wartość tą otrzymano przez pomnożenie liczby wyemitowanych akcji przez średnią wartość akcji w przeciągu roku.
ceny <- l.out$df.tickers$price.open
zwrot <- diff(log(ceny))
liczba_dni <- length(zwrot)
zmiennK <- sd(zwrot)*sqrt(liczba_dni)
wartKap<-83000303*mean(ceny)
r<-0.0625
T<-1
zobow<-9765000000
#Funkcja licząca błędy pomiędzy wartością kapitału oszacowaną oraz rzeczywistą wraz ze zmiennością kapitału
MertonModel<-function(x){
V0<-x[1]
zmiennV<-x[2]
d1<-(log(V0/zobow)+(r+zmiennV^2/2)*T)/(zmiennV*sqrt(T))
d2<-d1-zmiennV*sqrt(T)
A<-V0*pnorm(d1)-zobow*exp(-r*T)*pnorm(d2)-wartKap
B<-pnorm(d1)*zmiennV*V0-zmiennK*wartKap
return(A^2+B^2)
}
# Funkcja minimalizująca błędy
min_blad <- function(x){
wynik <- optim(c(V0=x[1],zmiennV=x[2]),MertonModel)
return(c(wynik$par[1], wynik$par[2]))
}
# Prawdopodobieństwo bankructwa
PD <- function(V, zob, r, zmienV, T){
d1<-(log(V/zob)+(r+zmienV^2/2)*T)/(zmienV*sqrt(T))
d2<-d1-zmienV*sqrt(T)
return(pnorm(-d2))
}Za pomocą funkcji min_błąd, oszacowano rynkową wartość aktywów wraz z jej zmiennością (V0 i zmiennV). Za pomocą funkcji optim, zminimalnizowano błędy funkcji “MertonModel” powiązane z teoretycznymi i zaobserwowanymi parametrami.
wyniki <- min_blad(c(11000000000,0.0))
V0<-wyniki[1]
zmiennV<- wyniki[2]
pd <- PD(V0, zobow, r, zmiennV, T)
datatable(as.data.frame(data.frame(data=c(V0, zmiennV, pd), row.names = c("Wartość aktywów", "Zmienność", "PD"))))